Slow Cooling of an Ising Ferromagnet 
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Abstract. A ferromagnetic Ising chain which is endowed with a single-spin- flip Glauber dynamics 
is investigated. For an arbitrary annealing protocol, we derive an exact integral equation for the 
domain wall density. This integral equation admits an asymptotic solution in the limit of extremely 
slow cooling. For instance, we extract an asymptotic of the density of domain walls at the end of 
the cooling procedure when the temperature vanishes. Slow annealing is usually studied using a 
Kibble-Zurek argument; in our setting, this argument leads to approximate predictions which are 
in good agreement with exact asymptotics. 
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I. INTRODUCTION 

The quenching procedure that is used in numerous 
studies of coarsening is usually based on the following 
protocol: (i) Start at a high initial temperature > T c , 
where spins are disordered; (ii) Instantaneously cool the 
system to a final temperature that is lower than the crit- 
ical, Tf < T c . The most popular choice of the initial 
temperature is Tj = oo. Cooling to the critical tem- 
perature, Tf = T c , is also often studied, especially for 
one-dimensional systems where T c = or in quantum 
quenches. After such instantaneous quench, ordered re- 
gions begin to grow and it takes a long time (an infinite 
time in the thermodynamic limit) to reach the equilib- 
rium. Understanding of the dynamics of phase-ordering 
kinetics, particularly the scaling regime that develops 
long after the quench, is the main goal 

Although the very fast cooling is natural in a num- 
ber of settings, there are numerous applications (ranging 
from the formation of glasses to the expansion of the 
Universe) where slow cooling is appropriate. A popu- 
lar gradual cooling scheme posits that the temperature 
vanishes linearly in time, 



T(t) = To 



1 - - 

T 



(1) 



so as time varies in the interval (0,r), the temperature 
falls from Ti — Tq to Tf — 0. To characterize the cooling 
scheme (TTJ) or other gradual cooling procedure, it is nat- 
ural to focus on the properties of the state at the end of 
the cooling procedure, particularly to find the difference 
between this final state [for the cooling scheme ([1]), it cor- 
responds to t = t] and a ground state. If the cooling rate 
1/r is large, the problem is not interesting as during time 
interval (0, r) the system has hardly evolved. Interesting 
behaviors arise if cooling is very slow, r > 1, so that the 
system has evolved over long time interval (0, r) and its 
state is close to a ground state. Then the problem is to 
describe the deviation from the ground state. 

The limit of slow cooling has been investigated in nu- 
merous papers. Earlier work has been mostly driven by 
applications to cosmology and proposals to study analogs 
of the cosmological phase transitions in the laboratory 



d H. Recently, similar ideas and methods have been 
applied to the dynamics of slow quantum annealing, see 
e.g. [U-@]. A surprisingly little work has been done in 
the more standard context of thermal phase transitions. 
Among a few exceptions is Ref. [8j that studies the Ising 
chain supplemented with Glauber's spin-flip dynamics [9( 
and Ref. which investigates high-dimensional sys- 
tems where T c > 0, so that the critical temperature is 
passed during the quenching procedure. 

Overall, studies of slow quenching tend to rely on un- 
controlled approximations. This is not so surprising since 
little can be done analytically in the high-dimensional 
setting. Some of the uncontrolled approximations, e.g. 
the so-called Kibble-Zurek mechanism 0, rely on 
sound physical ideas, yet we still want to better under- 
stand their possible limitations. Numerical (see e.g. [ID - 
HH) and laboratory (see e.g. fl4T4lrj |) experiments help in 
this regard, but large size simulations can be challenging 
(particularly for systems exhibiting quantum phase tran- 
sitions), while real system always have features beyond 
those that are captured by theoretical models. 

Unfortunately exact results haven't been obtained even 
for the Ising- Glauber spin chain Q. Here we show that 
this model is tractable for essentially arbitrary cooling 
procedure. More precisely, we establish an integral equa- 
tion for the density of domain walls. This equation is 
closed and linear. Although it is generally unsolvable, in 
the interesting slow cooling limit it is possible to deduce 
analytically an exact asymptotic for the final density of 
domain walls. 

The following section [TT] we write-down the governing 
equations for the two-spin correlation function, outline 
the problems in solving these equations, and present our 
chief results. We then provide a complete derivation for 
a special cooling scheme (section IIIII) . In section IIVI we 
analyze an arbitrary cooling scheme and derive various 
asymptotic results announced in section [TT] Section [V] 
is devoted to the comparison of exact results with ap- 
proximate predictions following from the Kibble-Zurek 
argument. Finally in section I VII we give a summary 
and briefly discuss how the cooling protocol can affect 
the probability for high-dimensional Ising ferromagnets 
to (eventually) fall into a ground state. 
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II. COOLING OF A FERROMAGNETIC ISING 
CHAIN: MAIN RESULTS 



For this cooling scheme the temperature decays from in- 
finity (at t = 0) to zero (at t = r) according to 



In one dimension, T c = and hence to assure coars- 
ening we must cool the Ising chain to zero temperature: 
Tf = 0. The zero-temperature evolution will continue 
as long as there are domain walls, so if we let the Ising 
chain to evolve during the time interval t < t < oo, it 
will eventually reach a ground state. Therefore the den- 
sity of domain walls at the end of the cooling procedure, 
t = r, gives the quantitative measure of the distance 
between the state at the end of cooling and the ground 
state. The deriving of this final density is our major goal. 

For concreteness, we always assume that the initial 
condition is spatially homogeneous. Then the chain will 
remain (on average) spatially homogeneous. Therefore 
the two-spin correlation function depends only on the 
separation between the spins, Gk(t) — (<jAt)<jj+k{t)) ■ 
The two-spin correlation function satisfies [3, [l7| 



dGh 
dt 



= -2G k + j(G k _ 1 + G k+1 ) 



for k > 1. Relation 



G (t) 



(2) 



(3) 



plays the role of a boundary condition. Note also that 
the density p(t) of domain walls can be expressed via 
Gi(t). Indeed, the bond + 1) hosts a domain wall if 
i(l — SjSj+i) = 1 and therefore 



-(1 - <Ti<Ti+l] 



1 - d 



(4) 



The parameter 7 in @ is given by 7 = tanh(2/T). [We 
set the ferromagnetic coupling strength J to unity; other- 
wise we would have had 7 = tanh(2 J/T).] Even if we are 
mostly interested by the density of domain walls, there 
is no closed equation for this quantity. Mathematically, 
we need to solve an infinite set of coupled linear ordi- 
nary differential equations ©. Glauber solved Eqs. @ 
in the situation when 7 is a constant, see 0, QjJ. For the 
gradual cooling scheme, however, we must deal with the 
time-dependent parameter "f(t). Fortunately, a lot can 
be done analytically. 

The cooling scheme ([1]) is characterized by a simple 
linear time evolution of temperature, yet the governing 
equations ([2]) depend on time via parameter 7 which has 
a more complicated time dependence: 



7 = tanh 



2 



1 - 



(5) 



As a warm-up, in the next section we shall investigate an 
alternative cooling scheme with a linear behavior of the 
parameter 7 that appears in the governing equations: 



lit) = t/r 



(G) 



T(t) 



tanh 1 (t/r) 



(7) 



The cooling schemes ((TJ and ([7]) are morally equivalent. 
The formulas for the latter scheme are a bit simpler and 
for that procedure we can start with Ti = 00 when initial 
conditions are particularly simple. 

Thus for the cooling scheme ([7]) we must solve 

^ L = -2G fe + -(G fe _ 1 +G fe+1 ), fc>l (8) 
at t 

on the time interval < t < t. The initial conditions are 

G k (t = 0)=6 k>0 (9) 

and the boundary condition is (J3]). The analysis of 
Eqs. ([8]) shows that in the case of slow cooling, t> 1, 
the final density of domain walls is 



r(|) 



(10) 



Remarkably, an Ising chain supplemented with essen- 
tially arbitrary cooling scheme is tractable. For concrete- 
ness, we shall assume that the temperature monotonously 
decreases from To to T(r) = 0; then 7 = 7(f) is a 
monotonously increasing function that varies from 70 to 
7(t) = 1. Actually we do not need to postulate that 
7(i) is monotonous. We consider the situation when 
7(i) varies from 70 = tanh(2/T ) to j(t) — 1; since 
7 = tanh(2/T), it must obey < 7 < 1, no other as- 
sumptions are required. The monotonicity of "/(t) re- 
flects that in typical cooling schemes the temperature is 
monotonously decreased, and having made such an as- 
sumption we are assured that 7 lies within the bounds 
7o < 7 < !• 

It turns out that in the limit of slow cooling, r ^> 1, 
only the behavior of j(t) near the end of the cooling 
procedure plays a role. For instance, in the case of the 
algebraic behavior, 



l-j(t)~A 1 



when t 



(11) 



the final density of domain walls depends only on two 
parameters A and a, namely 



P (t) ~ G(2r)- tt / 2 ( 1+tt ) 



with amplitude 



G 



HI' 


-2aY 


-1 


' A ' 


[ V2- 


-2a) 




1 + a 



(12) 



(13) 



An algebraic approach (fTTj) of 7(i) to unity is quite nat- 
ural. However, for the cooling scheme (UJ) the behavior of 
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j(t) near the end of the cooling procedure is exponential 
rather than algebraic 



1 - j(t) ~ 2exp 



4/T 
1 - t/n 



when t — y t 



This suggests to analyze a family of cooling laws with an 
asymptotic behavior of j(t) of the form 



1 — 7(t) ~ Bexp 



{l-t/rf 



when t — y t (14) 



with arbitrary positive parameters B,b,f3. In this sit- 
uation, the asymptotic behavior of the final density of 
domain walls is 



(15) 



This is a remarkably universal asymptotic: The leading 
t -1 / 2 algebraic factor does not depend on the parameters 
B,b, P and the logarithmic prefactor depends only on one 
parameter /3. 

We now turn to derivations of the announced results. 



III. SPECIAL COOLING SCHEME 

In this section we consider the special cooling proce- 
dure ([S])- First, we employ an exact analysis and then 
turn to the asymptotic behavior of the final density of 
domain walls. 



A. Exact Analysis 

To handle an infinite set of ordinary differential equa- 
tions ([8]) we recast it into a single differential equation 
using a generating function technique. Multiplying (J8J) 
by z k and summing over all k > 1 we find that the gen- 
erating function 



G(z,t)=J2G k (t). 



(16) 



fe>i 



satisfies 
dG 

at ' 



-2 + - (z + z- 1 ) 



G + -[z-G 1 (t)} (17) 



Solving (jTT]) subject to G(z,t = 0) = 0, we get 



G= ( ^'-[z~Gi(t')]exp 
Jo T 



2(t' — i) + 77 - 



where rj = (i 2 — t' 2 )/r. Recalling that the exponential 
is the generating function for the 

00 



factor exp 
modified Bessel functions 

y-l 



z + z 



we conclude that 



G k (t)= M-e^'-^lh-xW-Gxit')!^)] (18) 



Thus to determine G\ we need to solve an integral 
equation 

d(t)= / dt'-e 2 (*'-*)[7o(r?)-Gi(t'KiW] (19) 



Higher correlation functions G k with k > 2 are then ex- 
pressed via G\ and (fl8|) . 

The above results are exact. Now we turn to the most 
interesting case of slow cooling, t> 1, and compute the 
final density of domain walls. 



B. Final Density of Domain Walls 

Using (j4|), we re- write the integral equation (fl9|) in 
terms of the density of domain walls: 



1 - 



= 2p(t) + 2 f dt'-eW-V h(v)p{t') (20) 



We are mostly interested in the final density p(r). Hence 
we specialize (|2U|) to t = r. It is also convenient to change 
the integration variable, t' — > rj = (r 2 — t' 2 )/r. We arrive 
at the integral equation 



1 - - 



&?£(»?, t) [I {ri) - h{ri)] 



2p(r)+ f d V E(v,T)h (V) P (r^T^) (21) 



where we used the shorthand notation 
£(?7, r) = exp 2r ^-\/l — ry/ 



T — 1 



(22) 



The chief contribution to the integral on the right-hand 
side of (|2"H is gathered when 77 ~ yjr [see Eq. (|26]l below] . 

Therefore we can replace p (r^/l — rj/rj by p(r) and the 

exponential factor £(77, r) by 



£(»/,t) 



-J7— J) /4r 



(23) 



as it follows by expanding the term inside the square 
brackets in Eq. (|22|) . Thus the integral on the right-hand 
side of (I2T1) becomes 



p(r) / d V e-^ 2 ^I 1 (r ] ) 
Jo 



(24) 



We can further use the asymptotic formula for the Bessel 
function, 



when rj ^> 1, 



(25) 
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and replace the upper limit in (|24l) by infinity. Thus we 
obtain 

r°° r ( -) t 1 / 4 

p(r) d V e- r > I^{2^)-^= P {t) U > (26) 

Jo V47T 

This obviously dominates the second term 2p(r) on the 
right-hand side of (fT0|) . 

On the left-hand side of (|10l) both terms are compara- 
ble. The integral 



dr,E(r,,T) Mv) ~ h(jl)] 



is gathered in the region 77 = 0(1). Thus we can replace 
£(77, r) by e _T) and integrate up to infinity. Using 



d V e- r '[I a (v)-h(v)} = 1 



(27) 



we conclude that the left-hand side of (|2lj) is asymp- 
totically equal to 1 — 1/2 = 1/2. Equating this to the 
asymptotic expression (1261) for the right-hand side we ar- 
rive at the announced results (|10j) for the final density of 
domain walls. 



A. Algebraic Cooling 

For the family of cooling procedures (fTTj) we have 

n 



2 / dt"j(t") 

Aft' 
= 2r-f )-—- r 1-- 
1 + a \ t 



l+a 



or equivalently 



.4 



1 _^ = ± 

t 2t 1 + a 



1+0 



which leads to 



e 2 (*'— ) = exp 



-V 



l + a (2r) Q 



(31) 



Using pip and (l2"5t we find that the integral on the right- 
hand side of Eq. (|29[) is asymptotically 



drj 




exp 



A 77 



l+a 



l + a (2r) Q 



(32) 



Computing the integral in (|32p we find that the asymp- 
totic of the right-hand side of Eq. (121)1) is given by 



IV. GENERAL COOLING PROCEDURE 

For a cooling scheme characterized by an arbitrary 
function j(t), the solution is found using the same pro- 
cedure as in Sect. MI Al Instead of the integral equation 
1(15)1 we get 

Gx{t) = [ A' T (i')e 2(t ' _<) ftfa) - G^)h{ri)\ (28) 



RHS 



2 / 3 + 2a 
,t V 2 + 2a 



1 + a\ 2 < 1 +°) 



.4 



(2r) 2 ( 1 +°)p(r) 



is asymptotically 



The left-hand side of Eq. 

LHS = 1 - \ ^ dr, e"" [J (rf) - h {rf)\ = 1 

where in the last step we have used ((2"T)). Equating the 
RHS to the LHS we arrive at (fl2l-(IT3l. 



with 77 = 2 J t , dt" j(t"). The final density of domain walls 
satisfies 

= 2p(r)+ [ drje^'-^Ix^pit') (29) 



Note that i' that appears in p(i') and e 2 ^ r ) should be 
expressed via 77 by inverting 



7,(0 = 2/ dt" 7 (t") 



(30) 



Further, the upper limit in integrals in Eq. (121))) is given 

byC = 7 7 (t' = 0)=2/ r dt" 7 (i"). 

The dominant contributions to both integrals in ((2"9")) 
are gathered in the region where 1 - f'/r « 1. In this 
region the relation between t' and 77 simplifies. We now 
perform the calculations for two families of cooling pro- 
cedures introduced in Sect. [TTJ 



B. Exponential Cooling 

For the family of cooling procedures (|14j) we have 

Computing the integral and inverting we find (keeping 
two leading terms) 



1 



n 



t 2t bp 

Hence the exponential factor in the integrands in (|29l) is 



e 2 (*'- T > = exp 



- r] -vp {£y^*p{-b(£ 



and therefore the RHS of ((2"9")) is asymptotically 
drj 



P(r) 



12-KT] 



exp 



■- (-) 
bB \2t) 



77 exp 
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The integral has a simple asymptotic behavior. Indeed, 
changing 77 — > £ via 77 = 2t(£/ lnr) 1 /' 3 we recast the 
exponential factor in the integrand into 



exp 



2B 



lnr 



1+1/0 



,-!-!>/« 



(33) 



In the r — > 00 limit, the exponential factor (|33p is asymp- 
totically equal to 1 when £ < 6 and zero when £ > 6. 
Therefore the RHS is asymptotically 



,2r(fc/ lnr) 1 /' 3 



2p(r)l 



6 
lnr 



(34) 



The LHS of (|29|) is asymptotically 1/2 as before. Equat- 
ing it to (|34|) we arrive at the announced result (fl"5|) . 



V. APPROXIMATE ANALYSIS 

A Kibblc-Zurck (KZ) argument has been proposed 
[I, 0] to determine the final density of defects (domain 
walls in our case). The KZ argument leads to qualita- 
tively correct results, yet it involves an uncontrolled ap- 
proximation and cannot predict more subtle quantitative 
features, e.g. amplitudes in the scaling laws. In most ex- 
amples, it is impossible to estimate the quantitative error 
as we do not possess (asymptotically) exact results. The 
ferromagnetic Ising chain provides a rare exception and 
hence it allows us to probe the quality of an approxima- 
tion provided by the KZ argument. 

The KZ argument suggests to estimate p(r) in the fol- 
lowing way: 

1. During the initial regime the density of domain 
walls p(t) is assumed to be given by its equi- 
librium value at the corresponding temperature, 
p(t) — p e q[T(t)]. This assumption is natural since 
the cooling rate is very slow. 

2. During the final regime the density of domain walls 
p(t) is assumed to be constant, p(t) = p* where 
P* = A=q(r*) with T* = T(t*). This freezing also 
seems plausible as the density is very small and 
during the final regime the density should not de- 
crease much. 

3. The Ising model endowed with a Glauber dynam- 
ics is characterized by the typical time r e q(T) to 
relax to equilibrium. We now postulate that the 
separation time is determined from the criterion 



First, we need to find the equilibrium solution. In other 
words, we must solve the stationary version of Eqs. 
These equations admit an exponential equilibrium solu- 
tion: By inserting Gk = rj k into @ we find that the 



right-hand side vanishes when 2 = 7(77 + 77 1 ). From this 
relation, r\ = j^ 1 1 — y/l — , and therefore 



I-77 _ 1 
2 ~~ 2 



1 _ WW 

7 



(35) 



The relaxation time is to a degree the matter of agree- 
ment. Equation $2$ suggests that 



1 



' eq 



2[1 - 7] 



(36) 



as a reasonable definition. Combining (|36l) and the cri- 
terion r cq (T») = t — t„ we arrive at 



2[l-7*] 



t — U, 7* = j(U) 



(37) 



Let us begin with the special cooling procedure ([T}. 
Then 7* = i*/r and therefore (j37|) gives 1 — 7* = l/y/2r. 
Hence 



7* 



2- 3 /V- 1 /4 (38) 



Comparing (jlO)) and (1381) we see that the exact amplitude 
C = Tr 1 / 2 /r(l/4) = 0.48887 is smaller than the predic- 
tion from the KZ argument C K z = 2" 3/4 = 0.59460. 



A. Algebraic Cooling 

For the family of cooling schemes ([TO we recover the 
correct scaling law (fT2")l . but with a slightly erroneous 
amplitude 



Ckz-2- 1 / 2 A 1 /(2+2«) 



(39) 



Interestingly, the KZ argument correctly predicts even 
the dependence on A: The ratio of the exact and ap- 
proximate amplitudes is the function of a alone 



C _ 0F 


K 3 " 


-2a\~ 


-1 


1 


Ckz 2 


L V2- 


~2u) 




1 + a_ 



(40) 



This ratio approaches to 1 as a -> and to \Fk]\ = 
0.8862269 when a -+ 00. The largest deviation from 
unity is 0.7913 . . ., so the maximal quantitative mistake 
of the KZ argument is about 20%. 



B. Exponential Cooling 

For the family of cooling procedures (|14[) we find 

1 /lnr\^ . , 



6 



Note that the ratio of the asymptotically exact density 
given by (fTo) to the prediction (j4"Tj) of the KZ argument 
is equal to Wtt/4; it is a universal number as it does not 
depend on the parameters B, b, j3 of the cooling proce- 
dures HU). The same ratio characterizes algebraic cool- 
ing schemes (ITT1) in the a — > oo limit. This seems to 
be a general feature: If ^(t) approaches to unity faster 
than any power law, then the KZ argument quickly gives 
a prediction and the exact answer is just \Jir/A times 
smaller. 



VI. DISCUSSION 

This work was inspired by the desire to test the validity 
of the Kibble-Zurek argument for the spin chain supple- 
mented with non-conservative (Glauber) dynamics. The 
conclusion is that in this setting, the KZ argument is 
very good, namely it gives correct exponents and even 
the amplitudes are approximated reasonably well. The 
chief outcome of this work is actually unrelated to the KZ 
argument or any other approximation, namely we have 
shown that the Ising-Glauber chain remains solvable in 
the general case of an arbitrary (externally imposed) tem- 
perature drive T — T(t). The term 'solvable' here means 
that for the most interesting characteristic, the density 
of domain walls, we have found a closed linear inhomo- 



geneous integral equation. It does not seem feasible to 
express a solution in quadratures, but interesting asymp- 
totically exact results have been extracted. We have also 
given expressions of other two-spin correlations functions 
via the density of domain walls. 

What would happen if we shall continue beyond the 
point t = r, that is, we shall follow the zero-temperature 
evolution for t > r. In one dimension, domain walls 
will continue to diffuse and annihilate upon collisions; 
eventually the system will reach a ground state. In the 
case of instantaneous cooling, higher-dimensional behav- 
iors are different: In two dimensions, the Ising ferromag- 
net may not reach a ground state (it can get trapped in 
a metastable stripe state)' in three dimensions, a ground 
state is never reached [18j . It seems plausible that for 
slow cooling the system may not get trapped, or the 
probability of getting trapped decreases as t — > oo. This 
appears an interesting challenge and in two dimensions 
at least, one may even hope for theoretical progress fl9j . 
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